tau = 15e-6;%s 
gamma = 1.0/tau; %s^(-1)
w0 = 2*pi*30e3; %s^(-1)
wc = 2*pi*30e3; %s^(-1)

x = 1;
bx = 0.1;
by = 0.1;
%%
tList=0:0.05:10;
[mx, my] =  mxmy(bx, by, tList, x, w0/gamma, wc/gamma);
plot(tList, mx, tList, my)
%%
f0List = 10e3:1e3:150e3;
acList =zeros(size(f0List)); asList =zeros(size(f0List));
fcList =zeros(size(f0List)); fsList =zeros(size(f0List));
for k = 1:length(f0List)
    f0=f0List(k);
    [ac, as, fc, fs] = amp_phase(x, 2*pi*f0/gamma, wc/gamma, 1);
    acList(k) = ac; asList(k) = as; 
    fcList(k) = fc; fsList(k) = fs; 
end
subplot(2, 1, 1)
plot(f0List, unwrap(fcList)/pi*180, f0List, unwrap(fsList)/pi*180)
subplot(2, 1, 2)
plot(f0List, unwrap(fcList)/pi*180 - unwrap(fsList)/pi*180)